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Abstract 

We analytically and numerically analyze groundwater flow in a ho- 
h— 5 mogeneous soil described by the Richards equation, coupled to surface 

I water represented by a set of ordinary differential equations (ODE's) on 

I parts of the domain boundary, and with nonlinear outflow conditions of 

Signorini's type. The coupling of the partial differential equation (PDE) 
i i and the ODE's is given by nonlinear Robin boundary conditions. This 

article provides two major new contributions regarding these infiltration 
conditions. First, an existence result for the continuous coupled problem 
is established with the help of a regularization technique. Second, we an- 
i^h alyze and validate a solver-friendly discretization of the coupled problem 

based on an implicit-explicit time discretization and on finite elements in 
space. The discretized PDE leads to convex spatial minimization prob- 
i " i lems which can be solved efficiently by monotone multigrid. Numerical 

experiments are provided using the Dune numerics framework. 

> 
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Figure 1: Typical shapes of the coefficient functions fc(s) and p c (s) 



1 Introduction 

This article is concerned with existence of solutions, and efficient numerical 
approximation of unsaturated groundwater flow with surface water and nonlin- 
ear in- and outflow conditions. To start with, let us introduce the considered 
model. Let ft C M. d , d = 1, 2, 3 denote a bounded domain occupied by a homo- 
geneous soil with boundary d£l of class C 1 and outer normal v. Let Ei n ,E out 
and Sat C dQ denote three relatively open, pairwise disjoint d — 1-dimensional 
C^-manifolds, representing the infiltration, outflow and Neumann boundaries, 
such that dfl = £; n U E out U Sn and Sj n n S ut = 0- The unsaturated ground- 
water flow is modelled in the time interval [0, T] using Richards equation (see 
e.g. [§]) for the water saturation s : ft x [0, T] — > [0, 1] and the water pressure 
p : fi x [0, T] -> K 

nd t s - div (fc(s)^ 1 (Vp + e)) = / inQx[0,T]. (1) 

Here n denotes the porosity, [i the viscosity, k the permeability, and e the 
gravity vector. For simplicity we set n = fi = 1 and e = (0, 0, 1) T in this paper. 
Due to our assumption of a homogeneous soil, the permeability k depends only 
on the saturation s. In the Richards model, the air pressure in the pore space is 
assumed to be constant. We consider it normalized to p gas — and replace the 
water pressure in ([T]) by the capillary pressure function p c (s). Thus, the capillary 
pressure p c has to be seen as a multivalued graph with p c (s) — {Pc{s)} for 
s < 1 and p c {l) = [p c (l),oo). Furthermore, we assume the residual saturation, 
meaning the saturation level below which no flow of water occurs, to be zero, 
such that p c is defined on [0,1]. This is necessary for the definition of the 
infiltration boundary conditions. Additionally, we suppose p c (l) = 0, which is 
essential for the limiting process on the outflow boundary. Typical shapes of 
k(s) and p c {s) are depicted in Fig. [T| It can be seen that the coefficient functions 
degenerate in the sense that fc(0) = and that k and p c may have slopes that 
are unbounded in the neighbourhood of s = 1 and s — 0, respectively. These 
degeneracies of the coefficient functions are one of the major challenges when 
treating ^ analytically. 

We impose three different types of boundary conditions. On E^r we prescribe 
homogeneous Neumann boundary conditions. On the outflow boundary E Q ut we 
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choose the Signorini-type boundary conditions 



qv>0, p<0, (q-v)p = on£ out x[0,T] (2) 

with q = —k(s)(Vp + e) the Darcy velocity of the fluid. This models a hy- 
drophilic porous medium that is in contact with air (cf. |31|). The inequalities 
([2| express that (i) water may not enter the domain, (ii) the pressure cannot 
be positive as p a ; r has been normalized to 0, and (iii) water can exit only if 
p = 0. Concerning the analysis given below, a major problem is to give sense 
to the product of traces in the last equality of Q . On the infiltration bound- 
ary Sj n , w e only consider inflow of water into the soil and ponding of water on 
the surface (without runoff). We assume that there is an external water source 
r(x, t),x £ X; n . The water at x either infiltrates directly into the soil, or accu- 
mulates in x (ponding). Hence the state of the surface water can be described 
by the water table w : E; n x [0, T] —> HL At the absence of surface water runoff, 
w is modelled as an ensemble of ordinary differential equations 

d t w = q ■ v + r on S in x [0, T], (3) 

which is coupled to the subsurface boundary flux q (to be understood in a 
suitable weak sense). Note that although Si n is the only part of the boundary 
through which inflow can occur, we may also have outflow through £; n . A 
possible model for the flux q on Sj n is to consider it as a vertical leakage through 
a very small semipervious layer of thickness b and hydraulic conductivity Kh |24| , 
which can be modelled (cf. [9, Chapt. 6.4]) as 

qv = . (4 

c 

Here c = b/Kh is the resistance of the semipervious layer. Its inverse c" 1 is also 
known as leakage coefficient. Unfortunately, Q is only valid for w > 0, as in 
the case of w = and a negative pressure p, we would obtain an inflow into the 
ground without water present on the infiltration boundary. Following [21] we 
modify Q by a w-dependent factor to obtain a new law valid for all w 

p + +p_min{l,(f) }-w 

qv= ± , (5) 

c 

where p + — max{p, 0} and p_ = min{p, 0}. Now, w = in ( |5| i mplies q ■ v = 
if p < 0. The parameter a is a regularization threshold, with (|5| being the same 
as @ for w > a. As we will see later, solutions of this model fulfill w > 
automatically. 

Existence results for unsaturated flow in porous media go back to the pio- 
neering work of Alt, Luckhaus, and Visintin [2]. The basis of their approach - 
which we will also use in this article - is the so called Kirchhoff transformation. 
[3"7] picks up the methods used in [5] and proves the existence of solutions for ([lj 
by regularization techniques, with a stronger solution concept than in [2]. The 
vital improvement is that [37J allows dry regions, meaning regions where the sat- 
uration lies below the (positive) residual saturation, on the outflow boundary. 
Hence, one has to modify the boundary conditions, as the pressure is not defined 
in dry regions. In |37| Robin boundary conditions are used in the regularized 
problem and convergence to the limit on the outflow boundary is obtained via 
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defect measures. Without considering outflow conditions, [TS] also studies dry 
regions and proves the continuity of the free boundaries between the saturated 
and unsaturated regions and between the unsaturated and dry regions. In [33 
uniqueness of the solution is proved for the problem introduced in [2 a by adapt- 
ing the methods developed in |32j . In order to show L^-contraction and with 
this the uniqueness of the solution, doubling of variables techniques are used as 
introduced in [21]. With respect to the condition on the infiltration boundary, 
in j2H condition (|5| is used to couple the Richards equation with a hyperbolic 
PDE modelling ponding of water and surface runoff on the infiltration bound- 
ary. Under the restriction d t s € existence and uniqueness is shown for 
this coupled system. We remark that we do not need this assumption in our 
approach. 

Concerning the numerical treatment of the Richards equation, a rich litera- 
ture can be found. For an overview we refer to |10[ Sec. 2.2] and [T3] and the 
literature cited therein. The Richards equation has been discretized using finite 
volume methods |221 127] , mixed finite clement methods [3] [35] , finite elements 
[251 125] . and discontinuous Galerkin schemes [5J |5U]. In most cases, the re- 
sulting algebraic systems were solved using Newton's method, which, however, 
suffers from ill-conditioning problems due to the degeneracies in the parameter 
functions (Fig. [IJ. In contrast, Berninger et al. [13] discretize the Kirchhoff- 
transformed Richards equation in a way that the resulting spatial problems 
can be solved efficiently using a monotone multigrid method. Since this ap- 
proach is based on convexity rather than on smoothness, it is also well-suited 
for non-smooth Signorini-type boundary conditions. Numerical studies in 1 13] 
demonstrate the efficiency of the solver, as well as robustness for extreme soil 
parameters which correspond to parameter functions that degenerate into step 
functions. An extension of this approach to heterogeneous soil using domain 
decomposition techniques can be found in (12.. Concerning the numerical cou- 
pling of Richards equation with surface water, we mention |2U1 ET1 155] . where 
shallow water equations are used as the surface water model. Although different 
discretizations are applied, all these approaches enforce continuity of pressure 
and normal flux, i.e., mass conservation, across the interface. These coupling 
conditions have also been considered in an approach based on the Kirchhoff 
transformed subsurface flow in [5] and |14| . where an implicit time discretiza- 
tion of the coupled subsurface and surface water models has been solved with 
a heterogeneous domain decomposition method. In [14], however, the pressure 
continuity was replaced by a leakage condition such as (|4]) that represented clog- 
ging of a nearly impermeable river bed. Accordingly, a Robin-Neumann type 
iteration was applied in [Hj, whereas in [8| a Dirichlet-Neumann type itera- 
tion was used to solve the coupled system. In the present article, we choose an 
explicit time discretization for the surface water as in |10[ Sec. 4.3]. 

Let us now give an outline of the article. In the first part of this article 
we prove a new existence result for the Richards equation endowed with 
nonlinear outflow conditions of Signorini's type ([2| and coupled to surface water 
by nonlinear Robin conditions |3j [5]) . The proof is done without requiring the 
assumption dts £ L 1 (J7). To achieve this goal, in Section|2]we introduce a global 
pressure u with the help of the Kirchhoff transformation and reformulate the 
system accordingly. In Section[3] we state the main existence theorem and derive 
its proof by combining two central approaches for existence and uniqueness 
results for nonlinear, degenerated parabolic PDEs, namely the regularization 
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technique (cf. [32 [37]) and the method of ^-contraction (cf. [2H 133] ). 

In the second part of the article, starting with Section [4] we present our dis- 
cretization of the coupled nonlinear problem that is implicit-explicit in time and 
uses finite elements in space. Via variational inequalities, the implicit-explicit 
time discretization leads to spatial convex minimization problems which con- 
tain nonlinear outflow and nonlinear physical Robin boundary conditions but 
no Dirichlet conditions. We prove unique solvability of these problems if the 
relative permeability is non-degenerate (e.g., regularized) and the water table 
is non-zero. In Section [5] we give a numerical example combining the coupled 
surface-subsurface problem and Signorini outflow conditions. We observe that 
the proposed time discretization leads to a stable method for reasonable time 
step sizes. We also observe that the surface water height remains nonnega- 
tive as predicted by the theory, with the exception of certain singular points. 
Our numerical approach proves to be useful even beyond the confines of the 
assumptions needed for the existence result. 



2 Global pressure formulation of the fully cou- 
pled system and assumptions on the data 

For a better assessment of the unsaturated flow equation, we use Kirchoff's 
transformation to obtain an equivalaent formulation for a global pressure u. 
We then summarize the resulting fully coupled system with nonlinear boundary 
conditions and detail all assumptions on the data that we will need for our 
analysis. 

Let us define the Kirchoff transformation 

PcO) , _ 



*(«):= k(p- l (q))dq, *(») = < -> ' ' (6) 

U$(l),Oo) if 8= 1, 







which yields V<I>(s) = k(s)Wp c (s). We then define the global pressure function 
u{x, t) 



®{s{x,t)) if s(x,t) < 1, 

*(l) + fc(l)(p(i,t)-p c (l)) if s(x,t) = l. 



We thus have Vit = k(s)Vp for a sufficiently smooth solution (s,p) of (|T|). 
Therefore, we obtain as equivalent versions of the Richards equation 

nd t s = yT 1 div(Vu + k(s) ■ je) + f(s), u G <f>(s), 

(7) 

or nd t $~ X {u) = fi^ 1 div(Vu + fc($~ 1 (u)) • 7 e) + /($" 1 (w)). 

Conversely, the inverse Kirchhoff transformation : u i— > p can be defined as 
follows (cf. [T3]) 

-i/ \ \-h{u) if u < 0, , . . f 1 

P = K ~ iu):= \i^ else" ' H{U) = J W^W) ^ (8) 

With the global pressure formulation (|7| we are now prepared to state our fully 
coupled system. 
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2.1 Fully coupled system for unsaturated flow in porous 
media with infiltration and Signorini-type outflow con- 
ditions 

Let us define tt T = x [0, T], fi = n X {0} and E M = S 4 x [0, t] for f G [0, T] 
with i = in, out, TV. With given initial values so and wq, we get the following 
system of differential equations for the unknown global pressure u and surface 
water height w 





s), and dts = div(Vu + k(s) ■ e) + /(a) 


in 




-(Vu4 


■ fc(s) ■ e) • i/> 


on 


S ou t,T 




u < #(1) and k 2 (s)s - fc 2 (l)l < 


on 


S ou t,T 


o < - 


((Vu + fc(s) ■ e) • i/) ■ (fc 2 (s)s - fc 2 (l)l) 


on 


S ou t,T 


(Vm4 


- fc(s) ■ e) • i/ = 


on 


^N,Ti 


(Vu4 


■ k(s) ■ e) ■ v = g(u, w) 


on 






s = s 


in 


no, 




d t w = (r - g(u,w)) 


in 


5jin,T) 




W = Wo 


on 





U £ 



AT,T) 



(9) 



where 

. . ?i+ w h(u)ip(w) , ( (w 

ff( w ' u ') = -fc(Y)^ + - + " > and W = min|l,[- 

As the mappings s i-> k 2 (s)s and s H > p c (s) are both monotone functions, and 
as Pgas = and p c (l) = 0, the outflow conditions stated in ^ and the outflow 
conditions prescribed in (J9j) are formally equivalent . 

It is also possible to allow p c (l) < as it is the case, e.g., for Brooks-Corey 
parameter functions [16 . But then the outflow conditions in ([2| and the outflow 
conditions in (J9j) are no longer equivalent. We would have to deal either with 
the conditions in (pi) or with similar conditions formed by replacing p by u in 
^ as done, e.g, in [1U|. In such case it would be very difficult to give a meaning 
to the traces in the last equality in as it demands an a priori estimate for 
|| div(g)||i2(Q), which in turn necessitates an estimate for ||9ts||/,2(m. In order 
to solve this problem, one can either assume ||9ts||£2m) < C or 3>'(s) > c\ > 0. 
If one makes one of these assumptions, it is also possible to show an existence 
result for the system with original outflow conditions based on the results 
in [53]. Otherwise, it is only feasible to derive a weighted L 2 -estimate for div(q) 
[57] and the original outflow conditions have to be replaced by some surrogate 
ones as suggested in [37] and Q. 

2.2 Assumptions on the data 

In order to guarantee among others the well-definedness of the boundary con- 
ditions, the coefficient functions have to fulfill several assumptions. Additional 



assumptions are needed for the proof of our existence result (see Theorem 3.1 1. 

The relative permeability and capillary pressure functions are supposed to 
fulfill the following assumption. 
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soil type 


fc(l)7//x [cm/day] 


a [cm l ] 


I I-} 


Hygiene sandstone 


108.0 


0.0079 


10.4 


Touchet Silt Loam G.E.3 


303.0 


0.005 


7.09 


Silt Loam G.E.3 


4.96 


0.00423 


2.06 


Guelph Loam (drying) 


31.6 


0.0115 


2.03 


Guelph Loam (wetting) 




0.02 


2.76 


Beit Netofa Clay 


0.082 


0.00152 


1.17 



Table 1: Typical values for k(l)j/fi, a and I used for the van Genuchten 
parametrization (40] for five different soil types. 



Assumption 2.1. The permeability 

k 6 C 1 ([0, 1], [0, oo)) is monotonically non- decreasing 

and k(s) — for s = 0. The capillary pressure p c : (0, 1] — > {0, 1} R is a 
monotone graph given by a function 

p c G C 1 ((0, 1),R), monotonically increasing. 

We set p c (s) — {Pc(s)} for s <= (0,1) and p c (l) — [p c (l),oo). Furthermore, we 
suppose Pc(l) — (hydrophilic case). Finally, we assume that with a constant 
cq > the following conditions hold: 



d s Pc(s)>l/c Vs 6(0,1), 

\d s k(s)\ 2 <c k{s) VsG(0,l), 



k(s)\ Pc (s)\ + ^/k!7)d s p c ( S ) < c Vs € (0,1/2), 

(l-s)y/d sPc (s) <c Vs 6(1/2,1). 



(10) 



In order to guarantee that the condition on the infiltration boundary is well- 
defined, we additionally presume that the coefficient functions are chosen in 
such a way that 

u-eL°°(n T ) =^p-€L°°(n T ). (ii) 



Taking Assumption 2.1 for granted, we can prove it_ € L°°(fly). Together 
with (111 this yields the well-definedness of the conditions on the infiltration 
boundary. 

To justify the above assumptions from a physical viewpoint, we check if they 
can be fulfilled for the van Genuchten model |30] of the coefficient functions 

k(s) = S 1/2 • (1 - (1 - s 1 /" 1 )" 1 ) 2 , Pc{s) = --{S- 1 '" 1 - l) 1 /', 

v ' a 

with m = 1 - l/l. In Table [l] some typical values for fc(l)7//i — the hydraulic 
conductivity at full saturation — and two variable parameters a and / for five 
different soil types are listed. The values are extracted from [40_. Starting with 
the implication (11), we see that this holds true for / < 3. Concerning estimates 
( 10 1, we observe that they are all fulfilled except for the second one. This is due 
to the fact that the description of the relative permeability k(s) proposed by 
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van Genuchten is bounded for 8 = 1, but on the other hand its derivative blows 
up. However, since unbounded derivatives of k are hydrologically unrealistic, 
we can modify k(s) for s near 1 so that all estimates (10 1 are satisfied. We 
emphasize the importance of the third inequality in (10 1, which guarantees that 
k declines fast enough for small s to absorb the pressure p c , which goes to — oo 
in this case. 

Another parametrization of k{s) and p c (s) is given by Brooks and Corey [16], 
see also [T7] and [3D] . One can check that all estimates (fTob are satisfied for 
this parametrization, however, with regard to Assumption |2 . 1 1 we always have 
p c (l) < 0. Moreover, it is a characteristic property of the Brooks-Corey func- 
tions that the crucial implication (111 is never fulfilled. On the contrary, the in- 
verse Kirchhoff transformation (|8| is always ill-posed around the minimal global 
pressure u c = $(0) with k _1 (w c ) = — oo, cf. [TtJl Sec. 1.3]. With regard to the 
analysis in physical variables, we can consider regularizations of Brooks-Corey 
functions which exhibit the non-degenerate case k(s) > cq > 0, for which impli- 
cation (111 is satisfied, cf. [TPl Sec. 1.4.3] and Section [4] However, even though 
the numerics for the Richards equation is challenging for the degenerate Brooks- 
Corey functions, we use them in our numerical examples in Section [5] analogue 

to pi]. 



Concerning the source term, we assume that / : Q,? x [0, 1] — > K is bounded, 
continuously differentiable in all entries, and monotonically decreasing in s. 
Moreover, / shall satisfy f(x,t,l) < and f(x,t,0) > for all (x,t) G fix- 
Eventually, we suppose that s i-> f(x, t, s) is affine on (0, A) for small A. Further- 
more, we assume that r : £i n ,T —> K + is bounded and Lipschitz continuous. The 
initial conditions shall be given by functions so : fl ^ [0, 1] and wo : Si n — ► K + , 
with s G and w £ L°°(X: in ) and 

s(x, 0) = sq(x) a.e. in and w(x, 0) = Wo(x) a.e. on T,- m . 

We impose the following compatibility condition for sq: 

3p : Q -> E, p G Pc(sq) a.e. on {k(s ) > 0}, 

p €H 1 ' i {(l)nL oo (Sl), s | Ein >0. 



3 Existence result for the degenerate coupled sys- 
tem 

The goal of this section is to prove an existence result for system Q. Unfortu- 
nately, ([9]) exhibits some unpleasant properties. First, the system is degenerate 
in the sense that k(s) —> for s — > and that the derivative p' c (s) may be 
unbounded in the neighbourhood of s = and s = 1. As a result, the solution 
of the Richards equation lacks regularity, making the treatment of the boundary 
conditions challenging, as it is difficult to give sense to the appearing traces of 
the functions. Moreover, we deal with a coupled system due to the boundary 
conditions at the infiltration boundary. To achieve our existence result, we thus 
first regularize the coefficient functions k(s) and p c (s) with a small regularization 
parameter 5 in order to obtain a non-degenerate system. Afterwards, we fix S 
and decouple the system. Existence of solutions for the decoupled sub-problems 
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then follows by standard arguments. By applying Banach's fixed point theorem, 
we get a unique solution (u§, w$) for the regularized coupled system system, for 
fixed S. The crucial part when proving a priori estimates independent of S is 
the derivation of a maximum principle for us and ws- This can be achieved 
by considering classical solutions of a regularized parabolic system. Finally, we 
pass to the limit in S and get a solution for our original problem (J9j) . Our main 
result summarizes as follows. 

Theorem 3.1 (Existence of a solution pair (u,w) of system (|9|). Let the as- 
sumptions from Section^ be fulfilled and let the pair (us,ws) be a sequence 
of regularized solutions in L°°(n T ) C\ L 2 ((0, T), H 1 ' 2 ^)) x C°([0, T], i 1 (S m )) n 
L°°((0, T), -L 2 (£j„)). Then there exists a subsequence (us,w$) and a solution 
pair(u,w) e L 2 ((0,T),H 1 ' 2 (ft)) xC°([0, T], L 1 (Sj„)) such that 

u s ^u in L 2 ((0,T),fl" 1,2 (n)), 
ws^w in C^^TIL 1 ^)), 

and («, w ) are solutions of ^ in the distributional sense. Additionally, w > 
holds a.e. on S in /or all t 6 [0, T]. 

In the following subsections we will prove Theorem |3.1| as indicated above. 



3.1 Regularization of the coupled problem 

At first, we regularize the coefficient functions. 

Assumption 3.2 (Regularized coefficient functions). The regularized coefficient 
functions should fulfill the following conditions: 

k s e C 1 ([0,1], (0, oo)) and p s e C°([0, 1],R) piecewise C 1 

are monotonically increasing. For S — > we have ks \ k uniformly on [0,1], 
PS Pc uniformly on compact subsets o/(0, 1) and ks(0) = S 2 . We suppose that 
[J s ps([0, 1]) = K and, for simplicity, that p c (X/2) — p$(l/2) holds. Finally, ks 
and ps shall fulfill the inequalities of Assumption \2A\ 



The Assumptions |2.1| and |3.2| allow the definition of a regularized global 
pressure via the Kirchhoff transformation 

u s (x,t) = $ s (s) = / ksip^iq^dq, 



and guarantee that <&s is bounded from below and that $5 — > $ uniformly on 
compact subsets of [0, 1). We set u c — $(0) as the minimal global pressure. 



Example 3.3. Let Assumption 2.1 be fulfilled and suppose that c\s < k{s) < 
C2S 2 on (0, 1) for constants < C\ < 02- Then the following regularization fulfills 
Assumption \3.2\ 

f Pc(S) + ^ ViGl0.fl, 
k s (s) :=5 2 + k{s) Vse [0,1], p s (s) := { Pc {s) Vs e {6, 1 - 6], 



,Pc(l - «) + £ =^p= a Vse (1-5,1]. 
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Note that all inequalities (10 1, which have to be satisfied by the regularized 
coefficient functions, except for inequality ( |10[ )3, are needed for the a priori 
estimates. If ( 10 (3 is not fulfilled, we would still be able to prove convergence 
of the approximate solution pair to a limit (it, if), but this limit would not 
necessarily solve our problem. 

Using the regularized coefficient functions, we get the regularized Richards 
equation 

d t sg = div(Vit5 + kg(sg) ■ e) + fg, ug = $g(sg), 

with the pressure pg — ps(ss) and fg — f(x,t,sg(x,t)). In a straightforward 
manner we obtain the conditions on the Neumann and the infiltration boundary: 

(Vita + kg(sg) ■ e) ■ v = gs(ug, w$) on E in ,T 
(Vug + kg(ss) ■ e) • v = on Ejv,t- 

Inspired by |37| . we impose the following nonlinear Robin condition on the 
outflow boundary: 

- (Vit 5 + ks(ss) ■e)-v= ^k s {ss)(ps)+. (12) 

In this way, we keep the idea behind the original ouflow conditions, as water 
can flow out only if p > and if p > we have a very large outflow due to 
the scaling 1/5. On the other hand, boundary condition (12 1 is more handy 
than the original one, as we do not have to deal with a variational inequality. 
Rewriting ( |l2"| with = — (Vug + ks(ss) ■ e) one gets 

qa v = k 5 (s s ) , 



which, similar to our condition on the infiltration boundary, can be interpreted 
as a pressure difference between the pressure inside the porous medium and the 
pressure outside p gas = 0. The initial conditions are replaced by 



s s (x,0) = s s (x) 



Pj'ipoix)) if s > 1/2, 
s (x) if s < 1/2, 



for all i6fl, guaranteeing ^(sg) € if 1,2 (fi). 

We get the following regularized system for the unknowns ug and wg 



u s = $s(sg), d t s s = div(ViM + kg(sg) ■ e) + f(sg) in fi T , 

-(Vug + kg(sg) ■ e) • v = -ks(sg)(ps)+, on E out ,T, 



(Vug + kg(sg) ■ e) • v = on Ejv.t, 

(13) 

(Vit«5 + kg(sg) ■ e) • v = gg(u Sl wg) on £ in ,T, 

/ n s s . \p5 1 {Po{x)) if s Q (x) > |, 

sg{x,0) = s := i \ m fi , 

[s (x) if s (x) < |, 

d t wg = (r - g s (ug,wg)) in E in ,T, 

wg = w on S ini o, 
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where 

(u s ) + w s h s (u s )ip(w s ) 



gs{u$,w s ) 



ks(l)c c 





hs(u S ) = / - 1 dr, ip(w s ) = min <^ 1, ( — J 

(u s )- 

3.2 Existence of a unique solution for fixed S 

Let 6 be fixed. To prove the existence of a unique solution pair (v,s,w$) £ 
L 2 ((0,T),.ff 1 ' 2 (fi)) x C°([0,T],L 1 (S in )) of we take the complete metric 
space M = C°([0, T], with the metric \\wi — w 2 \\m = max <t<T / s . | 

W2(t) \dx, construct a mapping A : M — ► M and show that this mapping is 
/c-contractive. The fixed point of the mapping A is the solution of the ODE- 
subproblem. Inserting this solution as data into the subproblem of the Richards 
equation results in the solution pair of the whole system. For the sake of read- 
ability, we omit the index 5 in this subsection. At first, we define the two 
subproblems. 

Definition 3.4 (Subproblem of the Richards equation). Let Wj-\ € M H 
L°°((0, T), L 2 {Yii n )) be given for j = 1,2,.... We characterize Uj as the so- 
lution of the problem 



dt^iuj) - div(Vu j + /c($" 1 (u J )) • e) = in fi 

-(Vttj + ft($ -1 («j)) ■e)-v=- fc($ -1 (u,-)) (pj) + on I 

on Sjvt 



(14) 

{Vuj + k(®~ l {uj))-e)-i> = g(uj,Wj- 1 ) on £ in ,T, 

$~ (uj)(x, 0) = sq in Qq, 



with 



g(uj,Wj-i) 



o 

h( Uj )= J k{ ^ 1{r)) dr, ^(«, J _ 1 )=min{l,(^i) + 

(%■)- 

Definition 3.5 (Weak formulation of the subproblem of the Richards equation). 

FFe say £/ia£ Uj G £ 2 ((0, T), 1,2 (£1)) is a weaA; solution of problem ( |14[ ), if 
$~ 1 (u J ) € L 2 ((0,T),L 2 (ft)) and 

$" 1 (M J )(0,x)(^(0,a;)da;+ / ^(u^dttp 



(Vuj + fc(*" i K)) • e) • Vy> + / f($- L { Uj )) <p 

^(k(®- 1 (u j ))( Pj ) + )v+ j gfaw-Jip (15) 
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for all test functions if E L 2 ((0, T), iJ 1 ' 2 ^)) with d t (p E L 2 ((0, T), L 2 (f2)) and 
<P<?) = 0. 

Lemma 3.6 (Existence of a unique weak solution) . LettVj-i E Mr\L°°((0,T), L 2 (Y, in )) 
be given for j E N, and let the assumptions from Section^ and Subsection \3.1\ 
be fulfilled. Then there exists a unique weak solution Uj E L 2 ((0, T), H 1,2 (Q)) 



in the sense of Definition 3.5 



Proof. A proof can be found, e.g., in [30J. □ 

Definition 3.7 (ODE subproblem) . Let Wj-i E M n L oo ((0, T), L 2 (£ in )) and 
Uj E i°°(r2 T ) n X 2 ((0,T),i? 1 ' 2 (fi)) 6e given for j = 1,2,.... We characterize 
Wj as the solution of the problem 



with 

9(uj,Wj,Wj-i) 



d t Wj = r - g(uj,Wj,Wj^i)) in E in<T 
Wj = w on S in)0 , 

fc(l)c c c 



(16) 



o 

h(Uj)= I k{ ^- 1{T)) dr, ^K_ 1 ) = min|l,( U ^i) + }. 

Lemma 3.8 (Existence of a unique solution of the ODE subproblem). Let 

Wj-i E M n L°°((0,T),L 2 (£*„)) and u, E L°°(0 T ) n £ 2 ((0, T), H^ 2 {n)) be 
given for j E N, and let the assumptions from Section [1] and Subsection \3.1\ 
be fulfilled. Then there exists a unique solution Wj E M of the initial value 
problem (16). 



Proof. It can be easily proved that the right-hand side of the ODE (16) fulfills 
the assumptions of the Caratheodory's existence theorem for ODEs, which in 
turn yields the claim. □ 

We proceed with some estimates for the functions uj and Wj. 

Corollary 3.9 (A priori bounds for fixed 5). Under the assumptions from Sec- 
tion [J| and Subsection \3.1\ and with a constant C > 0, the following estimates 
for fixed S > and arbitrary index j hold true: 

1. s 4 €[0,1], a s eL*([0,T\,H l >*{Sl)), J Ut \d t s s \ 2 < C, 

2. u s E L°°(n T )nL 2 ([0,T],H l ' 2 (n)), J Qt \d t u s \ 2 < c, 

3. || div(Vwa + ks(s s ) ■ e)\\ L 2 [nT) < C, 

I ps(s s ) E L°°{Q T ), p' 5 (s 5 ) E L°°{n T ), 

5. h s (u 6 ) E L°°(St), 

6. <&s{s$) is Lipschitz continuous in s$. 

Proof. This is a simple conclusion from the estimates for arbitrary 5, which will 
be derived in Subsection 13.31 □ 
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Now we can define our mapping A : M — > M as a composition of two 
mappings B : M -» L 2 ((0,T), H 1 ' 2 ^)) and D : M x L 2 ((0,T), H 1,2 (f2)) -> M 
with Uj := Bwj-i and lUj := D(wj—i,Uj), where B and D are defined by the 
subproblems ( 14 1 and ( |16[ ). So A can be defined as u>j = D(wj-i, Bwj-i) =: 
Awj-i. We now prove the fc-contractivity of A. 

Proposition 3.10 (Estimate for the boundary conditions and the saturation). 

Let Wj-i , u>j -i £ M f] L°°((0, T), L 2 (E ira )) be two solutions of the initial value 
problem ( 16 1 anduj,iij £ L°°(f2x) H-L 2 ((0, T), H 1,2 (fl)) the associated solutions 
of (14 1 for j = 1,2, . ... Furthermore, let the assumptions from Sectionlw and 
Subsection 3.1 be fulfilled. Then for < t < T, Sj = <& (v,j), Sj — $ (ij) 
and a constant C > the following estimate holds true: 



\ Sj {;t)-S j (;t)\ + 



1 



t k(l)c 



(Uj)+ - (Uj),\ + 



\h(uj) - h(uj)\ 



(17) 



-■\Hsj){p+)(sj) - k(sj)(p + )(sj)\ < C\\wj-i - u>j-i\\ L i iSint) . 



The constant C can be specified as C = c 1 (l + \\h(uj)\\L°o(£ in t ) L^) where 
is the Lipschitz constant oftp. 

Proof. To prove this estimate, we use duality techniques, which is done for a 
similar case in detail, e.g., in [23] and for the case here in [35]. In short, we 
put all derivatives onto the test function ip and solve classically a dual problem 
in tp, which, loosely speaking, subtracts the boundary conditions appearing in 
the weak formulation, and adds them in absolute values, arranged with the right 
algebraic sign for the estimate. □ 

From the integral formulation of the ODE we get with straightforward mod- 
ifications 



S Jn 



Wj \ < 



1 



k(l)c 
C 



(Uj)+ - (Uj) \ + 



S in , t 
Wj-l - t&j-l| 



(18) 



ip(wj-i) \h(uj) - h(uj)\, 



with C having been defined in Proposition 3.10 This estimate (18 1 together 

1 



with (17) give 



s ln 



\Awj-i - Auij-il = / \uij -Wj\ < 



S ln 



S in ,, 



\vjj — Wj \ + C 



S in , I 



\ W j~l - W J-l\ 



where C = 2C. The application of Gronwall's lemma yields 



wj\ < Ce< 



u j-i 



|W 7 _i - W~_ X \ 



< TCe " max 

0<t<T 



S hl 



which results in 



\\Aw_ 



3-1 



Aw 



■j-l\\M 



< TCe~' 



\ w i-l -Wj-i\\ M - 



(19) 



13 



For T-C-e? = < 1, Ais obviously /c-contractive. Since wq = Wq, ( 19 1 provides 
the uniqueness of the solutions Wj and Uj, too. With Banach Fixed Point 



Theorem we get existence of a unique solution pair (115, w$) of the regularized 
system (13 1. We summarize this result in the following theorem. 



Theorem 3.11. For each 8 > there exists a unique solution pair ug £ 
L°°(n T ) n L 2 ((0,T),H^ 2 (n)) and w s £ M n L°°((0, T), L 2 (E m )) of the reg- 
ularized system fT3| ). Furthermore, we have ws > a.e. on for all t £ [0, T] 
and each 8 > 0. 

Proof. It remains to prove that ws > a.e. on £i n for all t £ [0, T]. We discretize 
the regularized system ( 13 1 explicitly in time, where we omit the index 8 again. 



By w n and u n we denote the time-discrete approximation of w(-,t n ) and t n ) 
respectively for a time step t n £ [0, T], n = 0, 1, . . . , JV, to := 0, fjy := T and 
a time step size r = t n+ \ — t n . We show by means of mathematical induction 
that w n ( x) > a.e on £j n for all n £ N. As w°(x) > by Assumptions on 
we prove now the induction step w n (x) > > a.e. 

the scheme 



2.2 



the data 

on £j n . To this end, we consider for fixed ieS 



n+1 



(x) 



r(x,t n ) 



u n (x) + w n (x) M" n W)-min{l,(^) + } 
fc(l)c c c 

(20) 



for the solution of the ODE ( 16 ) on Sj n .T- Please note that ug(x) is a solution of 
a non-degenerate elliptic problem and therefore regular enough to be evaluated 



on E; n . Let w n (x) > a.e on S; n be fulfilled. Rearranging (20 1 leads to 

u n (x)+\ (h(u n {x))^{w n {x)) 



w 



n+1 



Or) 



w n (x) 



r{x,t n ) 



fc(l)c 



Next we check which values r guarantee that w n+1 (x) 
modifications one gets that the estimate 



> 0, too. By simple 



r < min ■ 



for x £ Sj n a.e. (21) 



. g -r(x,t n )+ h{uHx)y 
has to hold true if the second and third entry in the braces are nonnegative- 



othcrwise these terms do not appear in (21 1. Now, if we let r go to zero and 
use the global estimate 



r < ess min ■ 

te[o,T] 



h(u(x,t)) 



for x £ Si n a.e., 



we obtain ws > for all 8 > 0. Note that r € C°(Ei„, T ) (2.2) and u £ 



L°°(T,i n ^T) imply that the introduced explicit Euler method converges pointwisc 
a.e. on E,„. □ 



We remark that estimate (21 1 can in principle also be used for numerical 



simulations in order to guarantee that the numerical solution remains positive. 



For each time t n , using an approximation of u(-,t n ) in (21 ), an upper bound Dt 
for t can be computed with pi) ensuring that the approximation of w(-, t n+ i) 
remains nonnegative for all r < Dt. However, numerical tests show that this 
bound is quite strict, see Section [5] 
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3.3 A priori bounds 

In this subsection we derive a priori bounds for the solution pair (u$,ws) of 



the regularized system ( 13 1 that are independent of 5. Let the assumptions of 



Section [2] and Subsection 13 . 1 1 be fulfilled. 

Lemma 3.12. There exists C < 00 independent of S > such that, for all 
5>0, 

I k sP ' s \d t s 5 \ 2 <C and f k 5 (s 5 ) \ div(Vu« + k s (s s ) • e)| 2 < C. 



Proof. The first estimate is obtained by testing the weak formulation (15) with 
dtu, using Assumption |3.2| and writing the integrands of the integrals on the 
boundary as total time derivatives. The second estimate follows directly from 
the first one. For a detailed exposition of the proof we refer to [38; and, in the 
case of Dirichlet boundary conditions on the infiltration boundary, to |37| . □ 

Theorem 3.13 (Maximum principle). There exists C < 00 independent of S 
such that, for all S > 

II (ug)+ ||i~(s <n , T ) + ||w<s|U°°(£i„, T ) < C, and \\us\\L™(n T ) < C. 

Furthermore, there exists p max < 00 independent of 5 such that p$(x,t) <p max 
for almost all (x,t) E fir and all S > 0. 

For the derivation of a maximum principle for the solution pair (us,ws), 
we seize ideas of [21] ■ We can obtain L°°-estimates independent of 5 in the 
case where the pair (ug,ws) can be approximated by smooth functions (u s s ,Wg) 
which are solutions of a regularized parabolic system. For the sake of simplicity 
we assume that the boundary of fl is of class C 2 . If dfl is not smooth enough 
one can proceed as in |24| and require that the elevation of the surface can be 
smoothly and periodically extended to the whole M 2 . Furthermore, we presume 



for the proof of Theorem [3T3] that k s E C 2 ([0, 1], (0, 00)) and p s E C 2 ([0, 1],R), 
otherwise an additional regularization step has to be performed. We consider 
the parabolic system 

u| - *,(«!), d t s E s = div(Vuf + fc 5 (s|) • e) + / £ ( S |) in (l T , 

-(Vu e 5 + kg(s e s ) ■ e) • v = ^k & {s%)i e {p\)xe, on E ou t,T, 

(VUj + kg(s e s ) ■ e) • v = on E n ,t, 

k 5 (l)c 

(22) 



(V«S + *«) ■ e) • * = -P^ + ^ + on £ in , r , 



s E s(x,0) = s E 0S in Q, 

d t w £ s - eAw £ s = r £ + - -± SJ in E in , T , 

Vwf -^=0 on <9E in x [0, T] 
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where v-i is the outer normal of 9E; n and the regularized functions fulfill the fol- 
lowing assumptions. f £ (s, •) and f £ (-,x 7 t) shall be in C 2 (VLt) and C 2 ([0, 1]) re- 
spectively, f £ -> /uniformly on [0, l]xfl T and ||/ e ||c°([o,i]xQ T ) < ||/||c°([o,i]xfi T ) 
r e shall be in C 2 (tt T ), ||r £ || L oo (aT) < ||r|| io = ( n T ) and \\r £ - r\\ H i,2 { n T ) -> 
as e — >• 0. Furthermore, we require that £ £ g C 2 (IR) and that £ e is nonneg- 
ative and nondecreasing such that £ e (pf) = {ps)+ if P<5 € Mr U (e,oo) and 
l£ £ (pf)l < 1 + £ 1 an d analog for ug. x e € C^°(E out ) is a cut-off function with 
X £ = 1 if diet (a, <9E out ) > £■ Additionally, h £ shall be in C* 2 (M), h £ > 0, h £ = 
if w| > 0, h £l < and h £ h uniformly on M. Moreover, we require ip £ to be 
in C 2 (R), to be nonnegative and that ip £ — > ip uniformly on R + . Finally, the 
regularized initial values Wq s and Sq s shall be in C 2 (E; n ) or C 3 (ft) respectively 
with \\w^ s - too,j||ii( S)te ) -> and ||s £ j(5 - s ,5 1 1 z, 2 (fi) ~> as e -> 0. They also 
have to fulfill suitable compatibilty conditions. 



Lemma 3.14 (Existence of a smooth solution pair (ug,Wg) of system (22 1). 



Suppose that the regularized functions in system (22 1 fulfill the just mentioned 
assumptions. Then there exists a solution pair (tt|, u>f ) with u\ G C (TIt), Vu| G 
C°(rV), W f € C°(E~^),A W f e C°(S~^),a tW | e C°(E~^). Furthermore, we 
have s £ s £ C°(n T ) and Vsf G C°(fJ T ). 

Proof. The regularized coefficient functions fulfill all necessary assumptions in 
the respective theorems in Section 5 of [30 , which yield the required result. □ 

For the solution (it|, w £ s ) we can now derive a maximum principle and obtain 
L°°-estimates independent of e and 5. 

Lemma 3.15 (Maximum principle for the approximating solution pair (u|, u>f )). 



For the solution pair (ug,Wg) of 22 we have the following inequalities 

l|w!IU°°(£in, T ) - C and IK^D+IU^CEin.T) ^ max { £ > c ill w IIU~(s lrl , T ) +c 2 }, 

where the constants C, Ci and c 2 do no£ depend on 5 and e. 

Proof. The proof follows ideas presented in jSJ. We omit the indices S and e 
in this proof. Let u s G C 2 (f2) be the solution of 

Au s + V/c($ _1 (u s )) • e + F = in O, 

u s =0 on E in , 

(Vu s + k(^~ l (u s )) • e) ■ v = on Eat, 

(Vu s + fc($- 1 (u s )) • e) ■ v = -&- x £{u s )x e on E out , 

where F = \\f(x,t, s)||c°(f2 T x[o.i])- Analogous to |24| one can show that u s 
is nonnegative. Therefore the term Vfc($ _1 (u s )) is unproblematic as ( &^ 1 (w.f) 
tends to zero as 5 — > for uf > 0. Thus we can bound Au s + Vfc($ _1 (u s )) • e 
and see that hence u s s as 5 — > on E out . A classical maximum principle 
yields that HusllcMn) — with C independent of 5 and e. 

We fix t G (0, T), define m t := max^ T )e T,- T u +( x ' T ) anc ^ se f U a '■= u s (x) + 
m t for x E il. With the same techniques as in [1] one can prove that 

u(x,t) < U s (x) V(x,r) G O^. (23) 
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We have that 

u{x(i),i) = U s (x(t)) for some (x(i),t) G x [0,t] (24) 

with m t — u + (x(t),t). We suppose that m t > 0. otherwise the proof is easy. 
Then Q and ([24) yield \7U s (x(i),i) ■ u < Vu(x(i),i) ■ v, which in turn leads 
to 

^K) < Mx®j)_ _ {WUs{x{i)) + k{ ^-i {Us{xm . e) .„< v M +C2 , (25 ) 
k(l)c c c 

where ci does not depend on e and 6 and 

v(t) := maxWM, W(t) := max w(x, t). 

r£[0,t] x£Z in 

Then v is nonnegative, nondecreasing and Lipschitz continuous on [0,T]. For 
W (t) = w(x(t),t) with x(t) in Ei n we have Aw(x(t),t) < 0. At the points where 
W(t) is differentiable, there holds 

W , M< W( - T K\M\ i ?(u(x(t),t)) f(m t ) W(r) 

w (r) " — T + l|r|l ^-> + — fe(TF^ ^ T(iF " "V + a (26) 

At first we assume that v(t) = W(t) for some < r < t. Since then v(s) = W(t) 
for all s E [i~,t], we have v' = on (r,f). Thus let u(i) = W(t) for some 
f G (0,T). Then (|25j> and Q lead to u'(t) < C and hence < C for any 
£ G [0,T] and a constant C not depending on 8 and s. Therefore we have 

IKIU=o (EtoiT) < C (27) 



and with ( 25 1 and ( 27 ) also 



ll(wI)+IU~(s in , T ) < max{e,ci|K|| L =o (Sin T) +c 2 }. 



□ 



If we now pass to the limit in e in (22 1 and obtain the following estimates. 

Lemma 3.16 (L°°-estimates for ug and ws). For the solution pair (ug,w$) 
of (13) we have the following inequalities 

IMIl~(O t ) + ||«>i||L°°(£ in ,T) ^ C > ( 28 ) 
ll(««)+IU 00 (Sin,T) ^ C lH W 5|U~(Si„,T) + C 2: (29) 

where the constants C, C\ and C2 do not depend on S. 

Proof. First let us assume that s ,s G n H 12 (n), w j G L°°(S in ) n 

ff 1,2 (Si n ) and that they fulfill suitable compatibility conditions. Approximating 
so, s and wo, (5 by the initial conditions of ( |22[ ) with Sg ^ — » so,s in i? 1,2 (f2) and 
w o 8 w o,s m ■ff 1 ' 2 (Sin), we obtain for the solution pair of p2| the 

following a priori estimates 

ll«ilU»(n T ) ^ i>« <Pmai on Q, T , (30) 

|Vu|||i2 (nr) + ||5 t uf || L 2 (f2T) + ||9t4llL 2 (n T ) 



|div(V«f + fcj(sf) • e) 



L 2 (f2 T ) 



't^5llL 2 (S in , T ) 



< C(<5), 



max J|3i«'I(-,*)||ii(i: in ) + max ||Vw|(-,t)|| L i (s . n) 
te[o,T] v *e[o,T] 



< c(5)(i + ||a tW K-,o)|| L1(Sin) + ||Vu,!(-,o)|| L1(Sm) ), 



(31) 



(32) 
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where all constants C and p max are independent of e and the constants in (30) 



do not depend on S either. Here Lemma 3 in [37J yields ( 30 I and the estimates 
in (31 ) can be obtained by a testing procedure. Applying the vanishing viscosity 



method (see for instance |4|) results in (32 1. We can extract subsequences 



with Vuf 



w s (- 



Vu s in L' 2 (fl T ) and w e s (-,t) -s- w s (-,t) in L 1 (E in ) 
t) in _L 1 (Ej n ) for t £ [0, T] as e — ► 0, where (ug, ws) is the 



us in L 2 (Vt T ) 
with dtWg(-, t) 

solution pair of ( |T3| ). To verify that the limit of Wg fulfills the ODE in ( 13 1, we 
test with a function in Qf(I]j n ) and apply the fundamental lemma of calculus of 
variations to obtain that dtWg— r§— (us)+/(ks(l)c)+ws/c+(hs(us)tp(ws))/c = 
a.e. on £ in for all t € [0,T]. 
that (s 0>s , 



for all t € [0,T]. Thus (|28j and {[29]) follow 
,w 0t s) are only in L°°(Q) or L°° 



'(Si. 



Let us now suppose 
respectively. We approximate 



(so,<s, Wo t s) by C^°-functions (Sq S ,Wq S ) such that Sq 5 — > so,5 m i 1 (fi) and 
Wq S — > u>o t s in i 1 (Ei n ) as 77 — >• and denote the corresponding solution pair 
of the system ( 13 1 (sl,wl). Estimates (17) and (19l imply that uP s — > ws in 
C ([0,T];L^S in )) and s n s (-,t) -t s s (;t) in L 1 ^) for a - e - * G [0,^] as r/ 0. 



Moreover, (|Tt|) renders /, 



sa(-,t) in L 1 (f2) for a.e. t g 



n 



hg(u s )\ 

0. Thus a priori estimates (28 1 and (29), which hold true for 



hold true for (us,uig). 



as 
s ) also 
□ 



77 



Remark 3.17. Note that the parabolic maximum principle yields that ss G 
[0,1] a.e. in tt T (cf. /577 for details). 

Lemma 3.18. There exists C < oo independent of 5 such that, for all S > 0, 



ks(ss)\S/s s \ 2 < C and u E L 2 ((0, T), H 1 ' 2 ^)). 



(33) 



Proof. To prove ( 33 ) , one first proves for a constant C, independent of 5, the 
estimate 



ks{ss)\Vps\ 2 X{s s >i/2} + ks(ss)d s ps(ss)\Vs s \ 2 X{s s <i/2} < C, 



where xd denotes the characteristic function of a set D. To this end, we test 
the weak formulation (15) with suitable test functions. For ss < 1/2 we choose 
(sg — 1/2)- + 1/2 and for ss > 1/2 we pick r] + (ps) = (ps — p) 



+ ■ 



T)-(ss) = (ss — 1/2)- + 1/2 and for ss > 1/2 we pick ij+(ps) 
where ps(0) < p < ps(l/2). Theorem 3.13 and Remark 3.17 guarantee the 
boundedness of r)_(ss) and r/ + (ps). For more details we refer to |371 138j . □ 



3.4 Passing to the limit in 5 



be fulfilled 



Theorem 3.19. Let the assumptions of Sectionffiand Subsection \3.1 
and let (us,ws) be the unique solution pair of the regularized system (13 1. Then 
for a subsequence 8 — > there holds 



ss 

Ug 
Ws 



s 
u 



in L°°(n T ), 

m L 2 ([0,T],^ 2 (fl)), 

m C ([Q,TlL x (V in )), 



(34) 
(35) 
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for appropriate limit functions s : Qt [0, 1] 
The limits satisfy 



id w : X 



in,T 



u e $(s) in Q T , d t s = div(Vu + k(s) ■ e) + f(s) inV'(Q T ), 

-(V« + k(s) ■ e) • v > on SoutT U Sat,t 

u<#(1) and fc 2 (s)s - k 2 (l)l < on S otti)T , 

< ((Vu + k(s) ■ e) ■ i/) • 2 (s)s - fc 2 (l)l) on E out T , 

(36) 

(Vu + k(s) ■ e) ■ v = on Ejv,t, 

(Vu + k(s) ■ e) ■ v = g(u, w) on E injT , 

s = s in Qq, 

d t w = (r - g(u, w)) in E m , T , 

w = w a on Ej n0 , 



where 



l(u,w) = - 



h(u) 



u + 
k(l)c 



w h(u)ip(w) 

C C 



dr. 



ip(w) = min < 1, ( — 



(") + 



The traces in ( 36 I3 — ( 36 1 7 exist in the sense of distributions. Furthermore, there 



holds w > a.e. on Ej n for all t G [0, T]. 

Proof. The assumptions in Section [2] and Subsection |3 . 1 1 and the a priori bounds 
of Subsection lOyield d34| , (p) , s : fi T [0, 1] and fs -± f in L 2 {n T ). Affinity 



of f(x,t,-) on (0, A) for small A and the co mpact ness of s$ on (A, 1] results 
in fs — >• / in L 2 (57t)- Using again Remark 3.17 the assumption k(0) = 
and the compactness of s$ on (A, 1] provides ks(sg) —> k(s) in L 2 (Q,t) and 
Vks(ss) — ^ V/c(s) in (L 2 (S1 T )) 3 for 5^0. Therefore, we can identify (|36|i as 
the distributional limit of (IT3|i. 



Compactness of the families sg and us- We outline the ideas of the proof 
and refer to |37| for more details. As is unbounded, we only aim at compact- 
ness for ug away from regions with maximal saturation. To prove compactness 
for the familiy us , we use a sequence e — > and an associated sequence of cut-off 
functions a e <E C°°(Cl, [0, 1]) with a s (x) = 1 for all x <G Q with dist(x, dQ) > e. 
If we consider r) e (t;) := (^ + e)_ for fixed e. it can be shown that the familiy 
a e ■ r) e (us) is compact in L 2 (VIt-s) and that for a subsequence 8 — > there holds 



<5->0 



r) e (us) — > Vei u ) in L (fir) for any e > 0. 



(37) 



Here, the compactness of the familiy a e -r] £ (us) in L 2 (f2T- £ ) can be deduced us- 
ing the Riesz characterization of compact sets. The compactness of the familiy ss 
away from regions with zero saturation follows similarly. 



For the proof of relation ( 36 1 2 , the boundary conditions on the Neumann and 



the outflow boundary and the initial conditions we refer to [37] . We only remark 
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that (36)2-4, (36l6 and ( 36 Is follow more or less straightforward, while the proof 



of (36 1 5 is much more challenging and uses defect measures and a compensated 



compactness argument . We now turn to the proof of ( 36 1 7 . The problem we have 



to cope with on the infiltration boundary is the fact that h(u) is not Lebesgue- 
integrable for S := {(x, t) £ fly | s(x, t) = 0} not being a null set in the sense of 
Lebesgue. Since we are not able to bound hs(us) on S independently of 5, hg(ug) 
is not necessarily Lebesgue-integrable on S either. To prove convergence on the 
infiltration boundary, we exploit that p-(x, t) — —h(u(x,t)) and that the water 
pressure p(x,t) equals the capillary pressure p c (s 
we know that (— p). 

Before analyzing the limiting process for S — > 0, we show that our boundary 
condition on Sj n 



Thanks to Assumption 2.1 
is bounded and therefore that h(u) is bounded, too. 



(Vug + kg(sg) 



(ug) 



wg h s (us)ip{ws) 



fci(l)c 



is well-defined in the sense of traces. Since wg is only defined on E; n , we do 



not have to discuss the terms soley depending on wg. Due to Theorem 3.12 c), 
the function (u$)+ has a well-defined trace on Si n ,T, which leaves us the term 



hg(ug) and, particularly, because of Theorem 3.12 e) 



h' s (us) 



(1 - sign(u 5 )) 



2 {kg ($J l (l/2(« a - \ug\)))) 



For simplicity, we restrict our analysis to the areas where us < 0. This is 
not a constraint, because h' s (ug) is bounded in the complement of these areas 
in Sly anyway. We arrive at h' 5 (ug) ~ 1/ (kg ($>J 1 (ug))) . Obviously, h' s (ug) is 
unbounded if sg — ^J 1 (ug) — > 0. As already discussed at the beginning of this 
paragraph, this would be a contradiction to the fact that (— p)+ is bounded due 
to assumption ( [TT| ). This proves that h' s (ug) is in fact bounded and, therefore, 
that hg(ug) has a well-defined trace. 

Since V hg(r) e (ug)) is bounded in (L 2 (f2y)) 3 , we get a s ubsequence, converg- 



ing weakly to a limit h. Using (37) and Assumption 
prove that h 



2.1 



one can than easily 
h(r/ E (u)) in L 2 (f2y). If we now 



Vh(j] e (u)) and hg(rj £ (us)) 
account for the fact that hg(ug) only depends on the negative part of ug, we get 
that Vhg(ug) — ^ Wh(u) and hg(ug) — > h(u) in (L 2 (Sly)) 3 and L 2 (f2y), respec- 
tively, for S 0. This results in hg(ug) — > h(u) in L 2 (Si ni y). The convergence 



us 



u in L (Sin.r) implies directly the convergence (ug)~ 



in L 2 (E in , T ) 



for S — > 0. Applying these results and Gronwall's Lemma, one gets u>s — > w 
in C°([0, T], L 1 (E in )), implying that w x (0) — Wq and w > a.e. on E in for all 

<e[0,T]. □ 



4 Discretization and convex minimization 

In this section we introduce an implicit-explicit time discretization of Sys- 
tem ([£]) , followed by a space discretization with finite elements for the Richards 
equation. We are bound to this particular discretization in the sense that it 
allows for spatial convex minimization problems in the porous medium that on 
the discrete level can be treated by monotone multigrid. We refer to [13] for 
a detailed analytical and numerical study of this approach. Here, we extend 
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the approach in [13J, replacing Dirichlet conditions by nonlinear physical Robin 
conditions. For the treatment of Robin conditions, in which the physical pres- 
sure appears without nonlinear deformation, see also |T2] . The ODE's for the 
surface water are treated explicitly. For another example of coupling Richards' 
equation with a single ODE representing the water height of a lake given by 
a compartment model, we refer to |10l Sec. 4.3]. See [5] and [TJ] for implicit 
treatments of the coupling of Richards' equation with the shallow water equa- 
tions, where in the latter article — as in this paper — the coupling is also given 
by a Robin transmission condition. 



4.1 Time discretization 

We introduce our special implicit-explicit time discretization of System (|9| , first 
in a strong and then in a weak formulation. Special attention will be given to 
the nonlinear Robin condition on S; n . The Signorini-type condition will lead to 
a variational inequality. 



4.1.1 Strong formulation 



Given s n in f2 and w n as well as r n+1 on Sj n for time steps t n and t n +i, 
respectively, with r = t n+ i — t n , we define (s" +1 ,w n+1 ) at time step t n +\ as 
follows. 



,71+1 



£ $(s" +i ), s" +i =rdiv(Vu" +i + k(s n ) • e) + r/(s n+i ) + s 



n, 



,n+l 



< and - (Vu n+1 + k(s n ) ■e)-u>0 
0={(Vu n+1 + k(s n )-e)-v)-u 



71+1 



(Vw n+1 + k(s n ) ■ e) ■ v = 
(V!i' i+1 + k(s n ) -e) u = g{u 



71+1 „.,77\ 



W 



71+1 



71 + 1 



g(u n+1 ,w n ))+w n 



on E Q 
on E Q 

on S^, 
on E in , 
in S,-„. 



(38) 



Here, boundary condition (38)3, (38)4 is the discretization of the original 
Signorini-type condition ^ which is equivalent to ([9]) 3 ([£]) 5 ■ Boundary condi- 
tion ( 38 1 6 needs special consideration. First, we rewrite g on Ei n x E; n (defined 
in (|9|)) by inserting the inverse Kirchhoff transformation kT 1 : u 1— > p given in 



Equations (|8| and obtain 

-g{u 



■ n+ \w n ) 



IV" 

c 



,71+1 



h(u n+1 ) min{l,(^) + } 



fc(l)c c 
= c- 1 K - 1 (u n+1 ) + + c" 1 ^ 1 ^ 1 )- • i/)(w n ). 



Therefore, with p n+1 := n 1 (u n+1 ), boundary condition (38 1 6 reads 



(k(s n+1 )v P n+1 + k{ s n ) - e )-u + c -y; +1 + 



-l p n_+l 



(Vu n+1 + k(s n ) ■ e) ■ v + c- 1 K - 1 (u n+1 ) + + c- 1 ^- 1 ^ 1 ). • tp{w n ) = c^w' 1 

on E in and is a nonlinear Robin condition both in p n+l and in u n+ . The 
special structure of it when formulated in u n+1 stems from the fact that the 
real functions 



(39) 
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Figure 2: Left: degenerate = k^ 1 with $ _1 (0) = -1.3245 • \pb\ > -oo by 
Brooks-Corey, (bold line, data from Section[5]) and with weight %j;(w n (£)) = 0-2 
for u_ (dotted line). Right: regularized (k^)s with 6 2 — 0.1 and $ -1 (0) = — oo, 
(bold line, construction by kg(s) '■= max(fc(s), 8 2 ) as in |10[ Sec. 1.4.3]) and with 
weight r 4 ! (w n (£)) = 0.2 for w_ (dotted line). 



which, for £ G Ej n almost everywhere, induce a superposition operator on Ei n , 
are increasing, cf. |llj . This observation is essential for our way to determine 
u n+1 that we present below. In fact, for w n (£) > a, the function c-kJ is nothing 
but the real function kT 1 that induces the inverse Kirchhoff transformation 
Otherwise, the factor ip{w n ((;)) = ( — — )+ < 1 just leads to a smaller weight 
for negative values of p, see Fig. [2] Finally, note that system (38 1 only contains 



an explicit coupling of ground and surface water in the sense that the solution 
u n+1 in the porous medium depends implicitly on the known value w n from the 
previous time step while w n+1 can be explicitly determined once u n+1 is known. 

4.1.2 Weak formulation 

In contrast to the considerations in Section [3]our weak formulation of the spatial 
problem p8| is based on the second line of ([7]), concentrating on u n+1 instead of 
s n+1 as the unknown on f2. Recall that by the Kirchhoff transformation (|6j, we 
have u > 4>(0). Depending on the kind of degeneracy of k(p~ 1 (q)) for q — > — oo, 
we might have $(0) = — oo or — oo < $(0) < as in the case of Brooks-Corey 
parameter functions (see (51 1 and ( |52| or |10l Sec. 1.3] and Fig. Ejj. Anyway, 



for a weak formulation we should expect u n+1 to be an element of the convex 
subset 

K ■= {v G ff^n) : v > $(0), W| Sout < 0} 

of the Sobolev space ff 1 (il). Together with the Signorini-type outflow condi- 
tion Q, which is also encoded in K,, these constraints lead to the following weak 



interpretation of (38)i-(38)6 in terms of a variational inequality for u n+1 G /C 
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(for details consult |10l Sec. 1.5.2/3, 3.4.5], where e = —e z ) 



($- V 1+1 ) - rfi®- 1 ^ 1 ))) ■ (v - u n+1 ) + t I K Uu n+1 ) ■ (v - u n+1 ) 



T I Vu 

In 



n+l 



V(v-u n+1 )- f <5>-\u n )-{v~u n+1 ) 
Jn 



Tl fc($- i (u"))e- V{v-u n+i )-T / (v-u n+i )>0 

In Jx in c 



in 



Vv £ K. 



(40) 



4.2 Convex theory 

In this section, we show that under reasonable assumptions, the variational 



inequality (40 I is equivalent to a uniquely solvable convex minimization problem. 

Let x £ Q, and £ £ S in , be primitives of - rf(x, t n+ i($ _1 (-))), 
x e fi, and r«|, £ € respectively, on [$(0),oo) n R. Since - r/$ _1 
and are increasing (cf. Assumption 2.2), ^ x and Vl^ are convex so that the 
functional 



Jn Js in * 



<f>- L {u n ) -v + t / fc($- 1 (M ,l ))e- Vu-r / u 



Vw e AC (41) 



is strictly convex. Due to the assumptions on k, p c and /, the functions 
x £ f2, are bounded and Lipschitz continuous. Note, however, that we have 
k|(u) — oo for u \ $(0) in case of > so that as least for $(0) > — oo, 

the transformation induced by k|, £ € Ej„, is ill-posed in general and the second 
integral in (40 1 is not well-defined as an L 1 -integral, cf. left graphics in Fig. [2] 
However, if we regularize k and p c according to Assumption 3.2 the regularized 
functions (k^)s, £ £ S; n , turn out to be Lipschitz continuous on R with Lipschitz 
constant max M (/-cp^(u) = ^min{l, ( w ff )+} and the corresponding integral 



in (40 1 is well-defined, cf. right graphics in Fig. [2] The following theorem states 
unique solvability of the variational inequality ( |40| for regularized functions 
( K l)s, £ £ Sin, if their inverses have uniformly bounded Lipschitz constants for 
£ £ Sj n or j equivalently, if the weighting factor tp(w n (£)) is uniformly bounded 
away from zero on a submanifold of E; n with positive Hausdorff measure. One 



can even prove well-posedness of (40 1, cf. |10l Prop. 2.4.11] 



Theorem 4.1. Let <F, 



Assumptions 2.1 or 3.2 



x G f2, be induced by coefficient functions according to 
Let £ £ £j n , be induced by regularized coefficient 
functions according to Assumptions 3.2 and, accordingly, let «|, £ £ Sj n , &e 
replaced by (lii )s in (40 I. Furthermore, let w n £ L 2 (£i„). Then the variational 
inequality (40) is equivalent to the convex minimization problem 



e K : F(u n+1 ) < F(v) Vv e JC. 



(42) 



If w n > holds on a nonempty submanifold of T, in , then (42) is uniquely solv- 
able. 
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Proof. Both the equivalence proof and the proof for the unique solvability is a 
combination of the convex theory on the time-discretized Richards equation with 
Signorini-type boundary conditions presented in [101 Sec. 2.3] and |13| and with 
Robin conditions which can be found in [TIH Sec. 3.4.1] and |12) . Here, Lipschitz 
continuity of the regularized (k^)s is needed. The additional dependency of ^ x 
on x £ ft and on £ G S; n can be easily included in the theory, taking into 
account f(-,t n+ i,s n+1 ) £ L°°(f2) and ij){w n ) £ L 2 (S; n ) which, in particular, 
leads to well-definedness of all considered integrals. 

The only issue we need to consider here is the lack of a Dirichlet boundary 
with positive Hausdorff measure which guarantees coercivity of the functional F. 
In our case, the latter is provided by the nonlinear Robin boundary condition 
with the help of Poincare's inequality |19l p. 127] in which gives 



1 

W\ 



L 2 (fl) 



<cn\\Vv\\ L2{n) Vv£H 1 (n) 



(43) 



for some constant cq > depending on £7. In detail, for the nonempty sub- 
manifold Y, w of E; n where w n > holds, i.e., w n > wq for a wq > 0, we prove 
that 

IHIe. := ||V«|U»(n) + IHIl»(e«) Vt-Gtf 1 ^) (44) 

defines a norm on H 1 ^) which is equivalent to || • ||jji(q). To see this, note that 
by the trace theorem [151 p. 1.61], there is a cs„ > such that the estimate 

IMIe„ <CEjM| ffl(n) yv&H\n) 

holds for the seminoma || • \\s w . Moreover, if = 0, then we have v = c w € K 

on D, by (43 1 and, furthermore, v = c w on T, w by the trace theorem, i.e., c w = 0. 
Therefore, || • is a norm on i/ 1 (f2) equivalent to || • ||H 1 (n)) cf. [41 [ p. 153]. 

W.l.o.g., we assume ^(O) = for all £ G E in and finish the proof by showing 
that there is a eg > such that 



c S \\vf 



L 2 (£„ 



< 



"'Si,, 



We prove that c s v 2 < * 5 (u) for all ( € E K . For all £ £ S in , * e (0) is the 
minimum of * ? since ^(O) = r(«|)j(0) = by ^ and ([39]). Therefore, 
it remains to show that 2cgv < t(k^)s(v) for v £ R and all £ G T, w or, 
equivalently, that t(k^)^ > 2c$ for £ G By ([8]) and (39 I we can choose 

c 5 <^mirx{l,»}. □ 



4.3 Finite element discretization and monotone multigrid 

In the following, we shortly indicate our space discretization and numerical solu- 



tion technique for (40 1. For more details we refer the reader to |101 Sec. 2.5, 3.4.5]. 



We discretize (40 1 or (41), respectively, by piecewise linear finite elements. For 
simplicity, we consider the two-dimensional case of a polygonal domain Vl C M 2 . 
Let Tj, j € No, be a conforming triangulation of f2. The set of all vertices of 
the triangles in Tj is denoted by Afj , the set of those vertices lying on Sj n and 
S out shall be called Afj n and Af° ut , respectively. The triangulation shall resolve 
the parts of the boundary corresponding to different boundary conditions. 
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Let Sj C 1 (17) be the finite element space of all continuous functions in 
ii^ 1 (SI) which are affine on each triangle t € 7j. To discretize the surface water 
level w we consider the trace grid 7~ ln of 7} on Ei n , and its dual grid 7~ m '* . On 

T™'* we define the finite element space <Sj n as the space of all functions that 
are constant on each element. This is a suitable space for w, as the evolution 
equation (J3j) does not involve space derivatives. Note that each element of 7~ ln '* 
corresponds uniquely to a vertex q of AQ n which we call element center. We 
will therefore sometimes use vertices of A/"j n to label elements in Tj n '* . The 

nodal basis function corresponding to a node q e A/} shall be called \q*\ and 
the (piecewise constant) basis function corresponding to an element q of 7~ ln '* 

will be called Aq n ' (i) . 

The finite dimensional analogue of K is the convex set 

Kj := {v e Sj : v(q) > $(0) Vq E JVj A v(q) < Vq € JV° ut } . 

We define the ^-interpolation Ig j : C(S1) — > <Sj by evaluation at the vertices, 
and the ^^-interpolation I S in : C(Si n ) — > 5™ by evaluation at the element 

centers of 7~ m '*. This allows to define the finite element discretization of F 



in (41 ) by 



r'K) 'U + r / fc($- 1 (u"))e- Vv-r / • u Vu G /Cj, (45) 



Si, 



c 



where we assume that a discrete surface water level w" € 5j n is available and 
also used in the definition of (k^)s and, equivalently, of The discrete finite 
element solution it™ S 5j is assumed to be known, too. Then, with the weights 

:= / and ft* := f \f^\ 

we can note the following. 

Theorem 4.2. Assume that the conditions of Theorem ^. 1\ are satisfied. Then 
the discrete variational inequality for <E JCj : 

]T [$- 1 (^ +1 ( g ))-r/(ci>- 1 (^+ 1 (g)))] («(«) -^ +1 (g))^ 

+ r ]T (^) 5 (^ +1 ( ( z))( w ( (7 )-^+ 1 ( (Z ))/ l ;" 



r 



Vu] +L • V(« - u" +i ) - / • (v - u" +i ) 

si in 



+ t I fc($- 1 «))e ■ V(« - < +1 ) - r / • (« - < +1 ) > V« e /C, 



(46) 
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is equivalent to the discrete convex minimization problem 



u^ +1 eK 3 : Fj{u™ +1 ) < Fj(v) Vu € JCj. (47) 

If w n (q) > holds for a q £ Af" 1 , then ( |47[ ) is uniquely solvable. 

The proof of this theorem for fCj C Sj is analogous to the one of its con- 



tinuous counterpart, Theorem 4.1 for JC C For this we refer again to 



[TUl Sec. 2.5] and [T5]. We only point out that the second integral in the def- 



inition (45 1 of Fj is the same as if S) n is replaced by the space of piecewise 



linear finite elements on 7" m . As for (40), well-posedness can also be established 



for (46 1. Under reasonable assumptions, one also obtains a convergence result 



for the discretization. 



Theorem 4.3. Let the conditions of Theorem J^.l be satisfied and Tj, j > ; 
be shape regular and hj — maxtg^. diamt — > for j — > oo. Let w" £ L 2 (Ei„), 
j > 0, be given with w" — > w n in L 2 (Ei„) for j — > oo. Furthermore, we assume 
C°°(fi)n/C to be dense inJC. We define s n+1 = $~ 1 (u" +1 ) and, with Kg obtained 

we define p n+1 = 



3.2 



by regularization of n according Assumption 
(analogously, s™ +1 o,nd p™ +1 are defined). Then, for j 
convergence 



P 



n+l 



-> p n+1 in H X (tl) and I Sj p] +1 -> p" 



In addition, we have 

n+l . n+l 



in H\n) and I Sj s n 3 



n+l 



n+l 



-> oo, we have the 
(48) 

m L 2 {Q). (49) 
in L 2 {fl) (50) 



if $ 1 : ($(0),oo) — > M is Lipschitz continuous for the first and Holder contin- 
uous for the second statement, respectively. 



The proof of (48) is an adaptation of the one given in fTt?l Thm. 2.5.9], 



where we need to use the norm defined in (44 1 instead of the 7? 1 -seminorm as 



an equivalent norm in the solution space, which is i/ 1 (SI) here. For the proofs 
of @ and d50l see [13]. 



In terms of finite elements, a straightforward discretization of (38 1 7 reads 



w ■ 



,n+l 



in S" 



Here, we can assume that u™ +1 has been obtained by (46 1 and that the iterates 



, w™ £ SJ 1 are given. 
Quadrature of the second integral in (|46j is given by 5^ -interpolation of 



$ 1 (u 1 j). In addition, with regard to the numerical treatment of the discretized 
problem (46), it is necessary to use an upwind technique for the gravitational 



(convective) term fc($ _1 (u™))e in the third integral. In the finite element con- 
text, this can be achieved by adding an artificial viscosity term to the discretized 
convection, cf. [101 Sec. 4.2] and [T3"] . 



The discrete convex obstacle problem arising from ( 46 1 is defined on a set JCj 
which is a product of intervals. To such problems, monotone multigrid methods 
can be applied |28| . In short, these methods first minimize Fj by successive 
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Figure 3: Left: Setting of the numerical example. Right: Coarse grids. The 
actual grids were created from these grids by four steps of uniform refinement. 



one-dimensional minimization along the A^ -directions for q e A/}. After appli- 
cation of this nonlinear Gauss-Seidel fine grid smoother, a suitable coarse grid 
correction given by constraint Newton linearization in the smooth regimes of 
the nonlinear system given by (46 1 accelerates convergence. Linear multigrid 
convergence rates are obtained asymptotically. For details concerning the appli- 
cation of this method to a Signorini-type and Robin problem like ( 46 1 we refer 
to [H Sec. 2.7, 3.4.5]. 

In the next section, we will apply monotone multigrid to an example with 
Brooks-Corey parameter functions which do neither satisfy the conditions for 
the existence theorem 3.1 nor for the theorems given in this section which require 
Kirchhoff transformations given by regularized parameter functions. However, 
even though the equivalence result of Theorem 4.1 does not hold if (k{)s is 
replaced by arising from the Brooks-Corey parametrization (since ( 40 1 is 
not well-defined), it turns out that (42 1 is still a well-posed convex minimiza- 
tion problem in this case. The same is tru e fo r the discretized problem (since 



^(O) < one can even prove Theorem 



4.2 



here) so that (47 1 is well-posed 



and can be solved by monotone multigrid methods |10[ Sec. 3.4.5] 



5 Numerical example 

We close the article with a numerical example involving the coupling of surface 
water to ground water, and a Signorini-type outflow condition. We conclude 
that our model and solution algorithm is well-suited for practical applications. 

The problem setting is depicted schematically in the left part of Fig. [3] As the 
ground water domain, we use a 10 m x 1 m rectangle whose lower left corner shall 
be situated in the origin of M 2 . Zero-flow boundary conditions are prescribed on 
the left, right, and bottom sides, with the exception of two 0.5 m stretches E out 
on the lower left and right, where we set the Signorini outflow conditions 
The entire upper side Sj n couples with the surface water model through the 
flux conditions (|3| and (|5|. The resistance c is set to 10 5 s (cf. E2]), and for the 
threshold parameter a we choose 0.02 m. We use sand as the porous medium, 
with soil parameters n = 0.437 (porosity), s m = 0.0458 (residual saturation), 
sju = 1 (maximal saturation), pi, = —712.2 Pa (bubbling pressure), A = 0.694 
(pore size distribution factor), and K — 6.66 • 10 _9 m 2 (absolute permeability), 
taken from [Ml Tables 5.3.2 and 5.5.5]. We set p, = 1.002 • 10~ 3 Pas for the 
dynamic viscosity of water. The parameters pt and A characterize the parameter 
functions p i-> s = Pc 1 (p) and s i-> k(s), according to Brooks-Corey (T6j 00] 
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and Burdine [T7] of the form 



s = p -i(p) = I s » + - s ™) fe) ' for ^ Pfc < (51) 
\i for p>Pb, 



and 



/ - \ 3 +x 
fc(s) = *T , se[s m ,s M }. (52) 

The domain is discretized with a uniform triangle grid on 161 x 17 vertices, 
giving a mesh size of h = v^/lGm. We construct this grid by four uniform 
refinement steps of a uniform 10 x 1 grid. As a by-product we obtain a grid 
hierarchy suitable for our multigrid solver. To discretize the surface water model, 
we use the dual grid 7~ m '* of the trace of the subsurface grid on Ej n . The grids 
are depicted schematically on the right of Fig. [3] For the implementation we 
used the Dune libraries |6j together with the domain decomposition module 
dune -grid-glue [7]. 

We start with the initial pressure — 2- 10 4 Pa in the domain, which, taking the 



material parameters and (51 1 into account, corresponds to the constant initial 
saturation 0.1401. The initial surface water height is set to a constant m, and 
we prescribe a constant rainfall of 30mm/h w 8.33 • 10~ 6 m/s on [5, 10] x {1}, 
and zero elsewhere. The evolution is simulated up to a final time T — 350 000 s, 
which equals roughly four days. For the time step size r we choose 100 s. This is 
well below the upper bound on the time step size implied by the CFL condition 

T h W = h nn 

00sup ae[(u] |fc'(s)| Kgg(3 + 2/A) ' 

cf. |101 p. 223], which in our setting evaluates to approximately r < h ■ 1.14 • 
10 3 s/m— 200 s. For each time step we solve the spatial problem (46 1 to machine 



precision by the monotone multigrid method described in ]13J. We omit the 
solver convergence rates and refer the interested reader to [13] for robustness and 
efficiency studies. While solving to machine precision usually means more work 
than necessary, it eliminates the effect of solver inexactness from the following 
observations. 

Fig. [4] shows a few snapshots from the evolution. Initially, surface water 
accumulates on the right, where there is rainfall, and infiltrates into the soil. 
A saturation front traverses the domain. When the right part of £1 is almost 
filled, a hydrostatic pressure distribution starts to build up there, and the water 
starts to infiltrate the left part of f2. Around time step 2429, the ground is fully 
saturated. Most of the water entering the domain leaves it through E out , while 
a small part exfiltrates on the left part of E; n , and leads to a rising surface water 
level there. At the last time step, the subsurface pressure field is contained in 
the interval [0 Pa, 2.53 • 10 4 Pa], and the surface water height is contained in 
[0.23 m, 1.89 m]. 



In a second step, we try to numerically assess the time step bound (21 1. The 
proof of Theorem [XTT] shows that a nonnegative surface water height w n at time 
step n implies a nonnegative surface water height w' n+1 at time step t n +i if the 



time step size r = <„+i — t n is bounded by (21 ), which states 



r < mm{c,0i(x),62(x)} for almost all x G Sj n , (53) 
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n = 200 



n = 400 



n = 600 



n = 800 






Figure 4: Snapshots of the numerical evolution at time steps t n — n-r — n - 100 s 
as indicated, with the top black line giving the surface water height, and isolincs 
of the pressure p at integral multiples of 2 kPa. 
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Figure 5: The three terms of the time step restriction, plotted as functions over 
the time step number; dotted line: used time step size r = 100s. 



where 



I oo else 



if this term is nonnegative, 



and 



— , , C °L , if this term is nonnegative, 
oo else. 



Note that by construction, the minimum in ( |53[ ) is always a finite nonnegative 
number. 

In Fig. [5]we have plotted c, mrri xe E in &i(x) and mm xe s in 02(x), with u n re- 
placed by u" for j = 4, as functions of the time step number n. We observe that 



our time step is below the bound (21 1 for almost all time steps. In accordance 
with the continuous theory we observe nonnegativity of the solution everywhere, 
with one notable exception. A surface water element with initial value w p — 0, 
but next to an element q with w q > will drop below zero if the subsurface 
pressure is negative there. This is a discretization effect. Unlike in the continu- 
ous case, the surface water levels w p and w q couple through a subsurface basis 
function. Since w q > and the subsurface pressure is negative, there will be 
a flux into the subsurface which will in turn suck water out of the element p, 
making w p negative. Through this effect, negative surface water heights down 
to —0.013 m where produced in our example. 



It is plausible to try to use the bound (|53|), with u n replaced by it™, as a 



time step control mechanism. To check whether this bound is sharp, we have 
recomputed our example with the time step sizes r = 50 s, 100 s, 200 s, 400 s, 
800 s, 1600 s, 3 200 s. Only at the last value of r we do observe instabilities 
presumably caused by violations of the CFL condition. As it turns out, the 
graphs in Fig. [5] are virtually independent of the time step size. Also, we do not 
see negative water levels (except for the discretization effect described above) 
for any of the time step sizes considered. Hence we conclude that the sufficient 
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condition ( 53 1 in the spatially continuous case is generally to tight when applied 
to fully discretized problems. 
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